* COPYRIGHT (c) 1977 AEA Technology
*######DATE 22 Feb 1993
C       Toolpack tool decs employed.
C       SAVE statements added.
C       MA28J reference removed.
C       ZERO made PARAMETER.
C
C  EAT 21/6/93 EXTERNAL statement put in for block data so will work on VAXs.
C
C
      SUBROUTINE MA28A(N,NZ,A,LICN,IRN,LIRN,ICN,U,IKEEP,IW,W,IFLAG)
C THIS SUBROUTINE PERFORMS THE LU FACTORIZATION OF A.
C
C THE PARAMETERS ARE AS FOLLOWS.....
C N     ORDER OF MATRIX  NOT ALTERED BY SUBROUTINE.
C NZ    NUMBER OF NON-ZEROS IN INPUT MATRIX  NOT ALTERED BY SUBROUTINE.
C A IS A  REAL ARRAY  LENGTH LICN.  HOLDS NON-ZEROS OF MATRIX ON ENTRY
C     AND NON-ZEROS OF FACTORS ON EXIT.  REORDERED BY MC20A/AD AND
C     MC23A/AD AND ALTERED BY MA30A/AD.
C LICN  INTEGER  LENGTH OF ARRAYS A AND ICN.  NOT ALTERED BY SUBROUTINE.
C IRN   INTEGER ARRAY OF LENGTH LIRN.  HOLDS ROW INDICES ON INPUT.
C     USED AS WORKSPACE BY MA30A/AD TO HOLD COLUMN ORIENTATION OF
C     MATRIX.
C LIRN  INTEGER  LENGTH OF ARRAY IRN. NOT ALTERED BY THE SUBROUTINE.
C ICN   INTEGER ARRAY OF LENGTH LICN.  HOLDS COLUMN INDICES ON ENTRY
C     AND COLUMN INDICES OF DECOMPOSED MATRIX ON EXIT. REORDERED BY
C     MC20A/AD AND MC23A/AD AND ALTERED BY MA30A/AD.
C U     REAL VARIABLE  SET BY USER TO CONTROL BIAS TOWARDS NUMERIC OR
C     SPARSITY PIVOTING.  U=1.0 GIVES PARTIAL PIVOTING WHILE U=0. DOES
C     NOT CHECK MULTIPLIERS AT ALL.  VALUES OF U GREATER THAN ONE ARE
C     TREATED AS ONE WHILE NEGATIVE VALUES ARE TREATED AS ZERO.  NOT
C     ALTERED BY SUBROUTINE.
C IKEEP  INTEGER ARRAY OF LENGTH 5*N  USED AS WORKSPACE BY MA28A/AD
C     (SEE LATER COMMENTS).  IT IS NOT REQUIRED TO BE SET ON ENTRY
C     AND, ON EXIT, IT CONTAINS INFORMATION ABOUT THE DECOMPOSITION.
C     IT SHOULD BE PRESERVED BETWEEN THIS CALL AND SUBSEQUENT CALLS
C     TO MA28B/BD OR MA28C/CD.
C     IKEEP(I,1),I=1,N  HOLDS THE TOTAL LENGTH OF THE PART OF ROW I
C     IN THE DIAGONAL BLOCK.
C     ROW IKEEP(I,2),I=1,N  OF THE INPUT MATRIX IS THE ITH ROW IN
C     PIVOT ORDER.
C     COLUMN IKEEP(I,3),I=1,N  OF THE INPUT MATRIX IS THE ITH COLUMN
C     IN PIVOT ORDER.
C     IKEEP(I,4),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN
C     THE L PART OF THE L/U DECOMPOSITION.
C     IKEEP(I,5),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN THE
C     OFF-DIAGONAL BLOCKS.  IF THERE IS ONLY ONE DIAGONAL BLOCK,
C     IKEEP(1,5) WILL BE SET TO -1.
C IW    INTEGER ARRAY OF LENGTH 8*N.  IF THE OPTION NSRCH.LE.N IS
C     USED, THEN THE LENGTH OF ARRAY IW CAN BE REDUCED TO 7*N.
C W REAL ARRAY  LENGTH N.  USED BY MC24A/AD BOTH AS WORKSPACE AND TO
C     RETURN GROWTH ESTIMATE IN W(1).  THE USE OF THIS ARRAY BY MA28A/AD
C     IS THUS OPTIONAL DEPENDING ON COMMON BLOCK LOGICAL VARIABLE GROW.
C IFLAG  INTEGER VARIABLE  USED AS ERROR FLAG BY ROUTINE.  A POSITIVE
C     OR ZERO VALUE ON EXIT INDICATES SUCCESS.  POSSIBLE NEGATIVE
C     VALUES ARE -1 THROUGH -14.
C
C
C COMMON AND PRIVATE VARIABLES.
C     COMMON BLOCK MA28F/FD IS USED MERELY
C     TO COMMUNICATE WITH COMMON BLOCK MA30F/FD  SO THAT THE USER
C     NEED NOT DECLARE THIS COMMON BLOCK IN HIS MAIN PROGRAM.
C THE COMMON BLOCK VARIABLES ARE AS FOLLOWS ...
C LP,MP  INTEGER  DEFAULT VALUE 6 (LINE PRINTER).  UNIT NUMBER
C     FOR ERROR MESSAGES AND DUPLICATE ELEMENT WARNING RESP.
C NLP,MLP  INTEGER  UNIT NUMBER FOR MESSAGES FROM MA30A/AD AND
C     MC23A/AD RESP.  SET BY MA28A/AD TO VALUE OF LP.
C LBLOCK  LOGICAL  DEFAULT VALUE TRUE.  IF TRUE MC23A/AD IS USED
C     TO FIRST PERMUTE THE MATRIX TO BLOCK LOWER TRIANGULAR FORM.
C GROW    LOGICAL  DEFAULT VALUE TRUE.  IF TRUE THEN AN ESTIMATE
C     OF THE INCREASE IN SIZE OF MATRIX ELEMENTS DURING L/U
C     DECOMPOSITION IS GIVEN BY MC24A/AD.
C EPS,RMIN,RESID  REAL/DOUBLE PRECISION VARIABLES NOT REFERENCED
C     BY MA28A/AD.
C IRNCP,ICNCP  INTEGER  SET TO NUMBER OF COMPRESSES ON ARRAYS IRN AND
C     ICN/A RESPECTIVELY.
C MINIRN,MINICN  INTEGER  MINIMUM LENGTH OF ARRAYS IRN AND ICN/A
C     RESPECTIVELY, FOR SUCCESS ON FUTURE RUNS.
C IRANK  INTEGER   ESTIMATED RANK OF MATRIX.
C MIRNCP,MICNCP,MIRANK,MIRN,MICN INTEGER VARIABLES.  USED TO
C     COMMUNICATE BETWEEN MA30F/FD AND MA28F/FD VALUES OF ABOVENAMED
C     VARIABLES WITH SOMEWHAT SIMILAR NAMES.
C ABORT1,ABORT2  LOGICAL VARIABLES WITH DEFAULT VALUE TRUE.  IF FALSE
C     THEN DECOMPOSITION WILL BE PERFORMED EVEN IF THE MATRIX IS
C     STRUCTURALLY OR NUMERICALLY SINGULAR RESPECTIVELY.
C ABORTA,ABORTB  LOGICAL VARIABLES USED TO COMMUNICATE VALUES OF
C     ABORT1 AND ABORT2 TO MA30A/AD.
C ABORT  LOGICAL  USED TO COMMUNICATE VALUE OF ABORT1 TO MC23A/AD.
C ABORT3  LOGICAL VARIABLE NOT REFERENCED BY MA28A/AD.
C IDISP   INTEGER ARRAY  LENGTH 2.  USED TO COMMUNICATE INFORMATION
C     ON DECOMPOSITION BETWEEN THIS CALL TO MA28A/AD AND SUBSEQUENT
C     CALLS TO MA28B/BD AND MA28C/CD.  ON EXIT, IDISP(1) AND
C     IDISP(2) INDICATE POSITION IN ARRAYS A AND ICN OF THE
C     FIRST AND LAST ELEMENTS IN THE L/U DECOMPOSITION OF THE
C     DIAGONAL BLOCKS, RESPECTIVELY.
C NUMNZ  INTEGER  STRUCTURAL RANK OF MATRIX.
C NUM    INTEGER  NUMBER OF DIAGONAL BLOCKS.
C LARGE  INTEGER  SIZE OF LARGEST DIAGONAL BLOCK.
C
C SEE BLOCK DATA FOR FURTHER COMMENTS ON COMMON BLOCK VARIABLES.
C SEE CODE FOR COMMENTS ON PRIVATE VARIABLES.
C
C
C SOME  INITIALIZATION AND TRANSFER OF INFORMATION BETWEEN
C     COMMON BLOCKS (SEE EARLIER COMMENTS).
C     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0)
C     ..
C     .. Scalar Arguments ..
      REAL U
      INTEGER IFLAG,LICN,LIRN,N,NZ
C     ..
C     .. Array Arguments ..
      REAL A(LICN),W(N)
      INTEGER ICN(LICN),IKEEP(N,5),IRN(LIRN),IW(N,8)
C     ..
C     .. Local Scalars ..
      REAL UPRIV
      INTEGER I,I1,IEND,II,J,J1,J2,JAY,JJ,KNUM,LENGTH,MOVE,NEWJ1,NEWPOS
C     ..
C     .. External Subroutines ..
      EXTERNAL MA30A,MC20A,MC22A,MC23A,MC24A
C     ..
C     .. Data block external statement
      EXTERNAL MA28J
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,AMAX1,MAX0
C     ..
C     .. Common blocks ..
      COMMON /MA28E/LP,MP,LBLOCK,GROW
      COMMON /MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK,
     +       ABORT1,ABORT2
      COMMON /MA28G/IDISP(2)
      COMMON /MA28H/TOL,THEMAX,BIG,DXMAX,ERRMAX,DRES,CGCE,NDROP,MAXIT,
     +       NOITER,NSRCH,ISTART,LBIG
      COMMON /MA30E/NLP,ABORTA,ABORTB,ABORT3
      COMMON /MA30F/MIRNCP,MICNCP,MIRANK,MIRN,MICN
      COMMON /MA30I/TOL1,BIG1,NDROP1,NSRCH1,LBIG1
      COMMON /MC23B/MLP,NUMNZ,NUM,LARGE,ABORT
      REAL BIG,BIG1,CGCE,DRES,DXMAX,EPS,ERRMAX,RESID,RMIN,THEMAX,TOL,
     +     TOL1
      INTEGER ICNCP,IRANK,IRNCP,ISTART,LARGE,LP,MAXIT,MICN,MICNCP,
     +        MINICN,MINIRN,MIRANK,MIRN,MIRNCP,MLP,MP,NDROP,NDROP1,NLP,
     +        NOITER,NSRCH,NSRCH1,NUM,NUMNZ
      LOGICAL ABORT,ABORT1,ABORT2,ABORT3,ABORTA,ABORTB,GROW,LBIG,LBIG1,
     +        LBLOCK
      INTEGER IDISP
C     ..
C     .. Save statement ..
      SAVE /MA28E/,/MA28F/,/MA28G/,/MA28H/,/MA30E/,/MA30F/,/MA30I/,
     +     /MC23B/
C     ..
C     .. Executable Statements ..
      IFLAG = 0
      ABORTA = ABORT1
      ABORTB = ABORT2
      ABORT = ABORT1
      MLP = LP
      NLP = LP
      TOL1 = TOL
      LBIG1 = LBIG
      NSRCH1 = NSRCH
C UPRIV PRIVATE COPY OF U IS USED IN CASE IT IS OUTSIDE
C     RANGE  ZERO TO ONE  AND  IS THUS ALTERED BY MA30A/AD.
      UPRIV = U
C SIMPLE DATA CHECK ON INPUT VARIABLES AND ARRAY DIMENSIONS.
      IF (N.GT.0) GO TO 10
      IFLAG = -8
      IF (LP.NE.0) WRITE (LP,FMT=99999) N
      GO TO 210

   10 IF (NZ.GT.0) GO TO 20
      IFLAG = -9
      IF (LP.NE.0) WRITE (LP,FMT=99998) NZ
      GO TO 210

   20 IF (LICN.GE.NZ) GO TO 30
      IFLAG = -10
      IF (LP.NE.0) WRITE (LP,FMT=99997) LICN
      GO TO 210

   30 IF (LIRN.GE.NZ) GO TO 40
      IFLAG = -11
      IF (LP.NE.0) WRITE (LP,FMT=99996) LIRN
      GO TO 210
C
C DATA CHECK TO SEE IF ALL INDICES LIE BETWEEN 1 AND N.
   40 DO 50 I = 1,NZ
        IF (IRN(I).GT.0 .AND. IRN(I).LE.N .AND. ICN(I).GT.0 .AND.
     +      ICN(I).LE.N) GO TO 50
        IF (IFLAG.EQ.0 .AND. LP.NE.0) WRITE (LP,FMT=99995)
        IFLAG = -12
        IF (LP.NE.0) WRITE (LP,FMT=99994) I,A(I),IRN(I),ICN(I)
   50 CONTINUE
      IF (IFLAG.LT.0) GO TO 220
C
C SORT MATRIX INTO ROW ORDER.
      CALL MC20A(N,NZ,A,ICN,IW,IRN,0)
C PART OF IKEEP IS USED HERE AS A WORK-ARRAY.  IKEEP(I,2) IS
C     THE LAST ROW TO HAVE A NON-ZERO IN COLUMN I.  IKEEP(I,3)
C     IS THE OFF-SET OF COLUMN I FROM THE START OF THE ROW.
      DO 60 I = 1,N
        IKEEP(I,2) = 0
        IKEEP(I,1) = 0
   60 CONTINUE
C
C CHECK FOR DUPLICATE ELEMENTS .. SUMMING ANY SUCH ENTRIES AND
C     PRINTING A WARNING MESSAGE ON UNIT MP.
C MOVE IS EQUAL TO THE NUMBER OF DUPLICATE ELEMENTS FOUND.
      MOVE = 0
C THE LOOP ALSO CALCULATES THE LARGEST ELEMENT IN THE MATRIX, THEMAX.
      THEMAX = ZERO
C J1 IS POSITION IN ARRAYS OF FIRST NON-ZERO IN ROW.
      J1 = IW(1,1)
      DO 130 I = 1,N
        IEND = NZ + 1
        IF (I.NE.N) IEND = IW(I+1,1)
        LENGTH = IEND - J1
        IF (LENGTH.EQ.0) GO TO 130
        J2 = IEND - 1
        NEWJ1 = J1 - MOVE
        DO 120 JJ = J1,J2
          J = ICN(JJ)
          THEMAX = AMAX1(THEMAX,ABS(A(JJ)))
          IF (IKEEP(J,2).EQ.I) GO TO 110
C FIRST TIME COLUMN HAS OCURRED IN CURRENT ROW.
          IKEEP(J,2) = I
          IKEEP(J,3) = JJ - MOVE - NEWJ1
          IF (MOVE.EQ.0) GO TO 120
C SHIFT NECESSARY BECAUSE OF  PREVIOUS DUPLICATE ELEMENT.
          NEWPOS = JJ - MOVE
          A(NEWPOS) = A(JJ)
          ICN(NEWPOS) = ICN(JJ)
          GO TO 120
C DUPLICATE ELEMENT.
  110     MOVE = MOVE + 1
          LENGTH = LENGTH - 1
          JAY = IKEEP(J,3) + NEWJ1
          IF (MP.NE.0) WRITE (MP,FMT=99993) I,J,A(JJ)
          A(JAY) = A(JAY) + A(JJ)
          THEMAX = AMAX1(THEMAX,ABS(A(JAY)))
  120   CONTINUE
        IKEEP(I,1) = LENGTH
        J1 = IEND
  130 CONTINUE
C
C KNUM IS ACTUAL NUMBER OF NON-ZEROS IN MATRIX WITH ANY MULTIPLE
C     ENTRIES COUNTED ONLY ONCE.
      KNUM = NZ - MOVE
      IF (.NOT.LBLOCK) GO TO 140
C
C PERFORM BLOCK TRIANGULARISATION.
      CALL MC23A(N,ICN,A,LICN,IKEEP,IDISP,IKEEP(1,2),IKEEP(1,3),
     +           IKEEP(1,5),IW(1,3),IW)
      IF (IDISP(1).GT.0) GO TO 170
      IFLAG = -7
      IF (IDISP(1).EQ.-1) IFLAG = -1
      IF (LP.NE.0) WRITE (LP,FMT=99992)
      GO TO 210
C
C BLOCK TRIANGULARIZATION NOT REQUESTED.
C MOVE STRUCTURE TO END OF DATA ARRAYS IN PREPARATION FOR
C     MA30A/AD.
C ALSO SET LENOFF(1) TO -1 AND SET PERMUTATION ARRAYS.
  140 DO 150 I = 1,KNUM
        II = KNUM - I + 1
        NEWPOS = LICN - I + 1
        ICN(NEWPOS) = ICN(II)
        A(NEWPOS) = A(II)
  150 CONTINUE
      IDISP(1) = 1
      IDISP(2) = LICN - KNUM + 1
      DO 160 I = 1,N
        IKEEP(I,2) = I
        IKEEP(I,3) = I
  160 CONTINUE
      IKEEP(1,5) = -1
  170 IF (LBIG) BIG1 = THEMAX
      IF (NSRCH.LE.N) GO TO 180
C
C PERFORM L/U DECOMOSITION ON DIAGONAL BLOCKS.
      CALL MA30A(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IDISP,IKEEP(1,2),
     +           IKEEP(1,3),IRN,LIRN,IW(1,2),IW(1,3),IW(1,4),IW(1,5),
     +           IW(1,6),IW(1,7),IW(1,8),IW,UPRIV,IFLAG)
      GO TO 190
C THIS CALL IF USED IF NSRCH HAS BEEN SET LESS THAN OR EQUAL N.
C     IN THIS CASE, TWO INTEGER WORK ARRAYS OF LENGTH CAN BE SAVED.
  180 CALL MA30A(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IDISP,IKEEP(1,2),
     +           IKEEP(1,3),IRN,LIRN,IW(1,2),IW(1,3),IW(1,4),IW(1,5),IW,
     +           IW,IW(1,6),IW,UPRIV,IFLAG)
C
C TRANSFER COMMON BLOCK INFORMATION.
  190 MINIRN = MAX0(MIRN,NZ)
      MINICN = MAX0(MICN,NZ)
      IRNCP = MIRNCP
      ICNCP = MICNCP
      IRANK = MIRANK
      NDROP = NDROP1
      IF (LBIG) BIG = BIG1
      IF (IFLAG.GE.0) GO TO 200
      IF (LP.NE.0) WRITE (LP,FMT=99991)
      GO TO 210
C
C REORDER OFF-DIAGONAL BLOCKS ACCORDING TO PIVOT PERMUTATION.
  200 I1 = IDISP(1) - 1
      IF (I1.NE.0) CALL MC22A(N,ICN,A,I1,IKEEP(1,5),IKEEP(1,2),
     +                        IKEEP(1,3),IW,IRN)
      I1 = IDISP(1)
      IEND = LICN - I1 + 1
C
C OPTIONALLY CALCULATE ELEMENT GROWTH ESTIMATE.
      IF (GROW) CALL MC24A(N,ICN,A(I1),IEND,IKEEP,IKEEP(1,4),W)
C INCREMENT GROWTH ESTIMATE BY ORIGINAL MAXIMUM ELEMENT.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (GROW .AND. N.GT.1) W(2) = THEMAX
C SET FLAG IF THE ONLY ERROR IS DUE TO DUPLICATE ELEMENTS.
      IF (IFLAG.GE.0 .AND. MOVE.NE.0) IFLAG = -14
      GO TO 220

  210 CONTINUE
  220 RETURN

99999 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'N OUT OF RANGE = ',I10)
99998 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'NZ NON POSITIVE = ',I10)
99997 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'LICN TOO SMALL = ',I10)
99996 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'LIRN TOO SMALL = ',I10)
99995 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE INDICES FOUND OUT ',
     +       'OF RANGE')
99994 FORMAT (1X,I6,'TH ELEMENT WITH VALUE ',1P,E12.4,/,20X,
     +       ' IS OUT OF RANGE WITH INDICES ',I8,' ,',I8)
99993 FORMAT (' DUPLICATE ELEMENT IN POSITION ',I8,',',I8,
     +       ' WITH VALUE ',1P,E12.4)
99992 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'ERROR RETURN FROM MC23A/AD')
99991 FORMAT (' ERROR RETURN FROM MA28A/AD BECAUSE ',
     +       'ERROR RETURN FROM MA30A/AD')

      END
      SUBROUTINE MA28B(N,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,IFLAG)
C THIS SUBROUTINE FACTORIZES A MATRIX OF A SIMILAR SPARSITY
C     PATTERN TO THAT PREVIOUSLY FACTORIZED BY MA28A/AD.
C THE PARAMETERS ARE AS FOLLOWS ...
C N      INTEGER  ORDER OF MATRIX  NOT ALTERED BY SUBROUTINE.
C NZ     INTEGER  NUMBER OF NON-ZEROS IN INPUT MATRIX  NOT ALTERED
C     BY SUBROUTINE.
C A      REAL/DOUBLE PRECISION ARRAY  LENGTH LICN.  HOLDS NON-ZEROS OF
C     MATRIX ON ENTRY AND NON-ZEROS OF FACTORS ON EXIT.  REORDERED BY
C     MA28D/DD AND ALTERED BY SUBROUTINE MA30B/BD.
C LICN   INTEGER  LENGTH OF ARRAYS A AND ICN.  NOT ALTERED BY
C     SUBROUTINE.
C IVECT,JVECT  INTEGER ARRAYS OF LENGTH NZ.  HOLD ROW AND COLUMN
C     INDICES OF NON-ZEROS RESPECTIVELY.  NOT ALTERED BY SUBROUTINE.
C ICN    INTEGER ARRAY OF LENGTH LICN.  SAME ARRAY AS OUTPUT FROM
C     MA28A/AD.  UNCHANGED BY MA28B/BD.
C IKEEP  INTEGER ARRAY OF LENGTH 5*N.  SAME ARRAY AS OUTPUT FROM
C     MA28A/AD.  UNCHANGED BY MA28B/BD.
C IW     INTEGER ARRAY  LENGTH 5*N.  USED AS WORKSPACE BY MA28D/DD AND
C     MA30B/BD.
C W      REAL/DOUBLE PRECISION ARRAY  LENGTH N.  USED AS WORKSPACE
C     BY MA28D/DD,MA30B/BD AND (OPTIONALLY) MC24A/AD.
C IFLAG  INTEGER  USED AS ERROR FLAG WITH POSITIVE OR ZERO VALUE
C     INDICATING SUCCESS.
C
C
C PRIVATE AND COMMON VARIABLES.
C UNLESS OTHERWISE STATED COMMON BLOCK VARIABLES ARE AS IN MA28A/AD.
C     THOSE VARIABLES REFERENCED BY MA28B/BD ARE MENTIONED BELOW.
C LP,MP  INTEGERS  USED AS IN MA28A/AD AS UNIT NUMBER FOR ERROR AND
C     WARNING MESSAGES, RESPECTIVELY.
C NLP    INTEGER VARIABLE USED TO GIVE VALUE OF LP TO MA30E/ED.
C EPS    REAL/DOUBLE PRECISION  MA30B/BD WILL OUTPUT A POSITIVE VALUE
C     FOR IFLAG IF ANY MODULUS OF THE RATIO OF PIVOT ELEMENT TO THE
C     LARGEST ELEMENT IN ITS ROW (U PART ONLY) IS LESS THAN EPS (UNLESS
C     EPS IS GREATER THAN 1.0 WHEN NO ACTION TAKES PLACE).
C RMIN   REAL/DOUBLE PRECISION  VARIABLE EQUAL TO THE VALUE OF THIS
C     MINIMUM RATIO IN CASES WHERE EPS IS LESS THAN OR EQUAL TO 1.0.
C MEPS,MRMIN  REAL/DOUBLE PRECISION VARIABLES USED BY THE SUBROUTINE
C     TO COMMUNICATE BETWEEN COMMON BLOCKS MA28F/FD AND MA30G/GD.
C IDISP  INTEGER ARRAY  LENGTH 2  THE SAME AS THAT USED BY MA28A/AD.
C     IT IS UNCHANGED BY MA28B/BD.
C
C SEE BLOCK DATA OR MA28A/AD FOR FURTHER COMMENTS ON VARIABLES
C     IN COMMON.
C SEE CODE FOR COMMENTS ON PRIVATE VARIABLES.
C
C
C
C CHECK TO SEE IF ELEMENTS WERE DROPPED IN PREVIOUS MA28A/AD CALL.
C     .. Scalar Arguments ..
      INTEGER IFLAG,LICN,N,NZ
C     ..
C     .. Array Arguments ..
      REAL A(LICN),W(N)
      INTEGER ICN(LICN),IKEEP(N,5),IVECT(NZ),IW(N,5),JVECT(NZ)
C     ..
C     .. Local Scalars ..
      INTEGER I1,IDUP,IEND
C     ..
C     .. External Subroutines ..
      EXTERNAL MA28D,MA30B,MC24A
C     ..
C     .. Common blocks ..
      COMMON /MA28E/MP,LP,LBLOCK,GROW
      COMMON /MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK,
     +       ABORT1,ABORT2
      COMMON /MA28G/IDISP(2)
      COMMON /MA28H/TOL,THEMAX,BIG,DXMAX,ERRMAX,DRES,CGCE,NDROP,MAXIT,
     +       NOITER,NSRCH,ISTART,LBIG
      COMMON /MA30E/NLP,ABORTA,ABORTB,ABORT3
      COMMON /MA30G/MEPS,MRMIN
      COMMON /MA30I/TOL1,BIG1,NDROP1,NSRCH1,LBIG1
      REAL BIG,BIG1,CGCE,DRES,DXMAX,EPS,ERRMAX,MEPS,MRMIN,RESID,RMIN,
     +     THEMAX,TOL,TOL1
      INTEGER ICNCP,IRANK,IRNCP,ISTART,LP,MAXIT,MINICN,MINIRN,MP,NDROP,
     +        NDROP1,NLP,NOITER,NSRCH,NSRCH1
      LOGICAL ABORT1,ABORT2,ABORT3,ABORTA,ABORTB,GROW,LBIG,LBIG1,LBLOCK
      INTEGER IDISP
C     ..
C     .. Save statement ..
      SAVE /MA28E/,/MA28F/,/MA28G/,/MA28H/,/MA30E/,/MA30G/,/MA30I/
C     ..
C     .. Executable Statements ..
      IF (NDROP.EQ.0) GO TO 10
      IFLAG = -15
      WRITE (6,FMT=99999) IFLAG,NDROP
      GO TO 70

   10 IFLAG = 0
      MEPS = EPS
      NLP = LP
C SIMPLE DATA CHECK ON VARIABLES.
      IF (N.GT.0) GO TO 20
      IFLAG = -11
      IF (LP.NE.0) WRITE (LP,FMT=99998) N
      GO TO 60

   20 IF (NZ.GT.0) GO TO 30
      IFLAG = -10
      IF (LP.NE.0) WRITE (LP,FMT=99997) NZ
      GO TO 60

   30 IF (LICN.GE.NZ) GO TO 40
      IFLAG = -9
      IF (LP.NE.0) WRITE (LP,FMT=99996) LICN
      GO TO 60
C
   40 CALL MA28D(N,A,LICN,IVECT,JVECT,NZ,ICN,IKEEP,IKEEP(1,4),
     +           IKEEP(1,5),IKEEP(1,2),IKEEP(1,3),IW(1,3),IW,W(1),IFLAG)
C THEMAX IS LARGEST ELEMENT IN MATRIX.
      THEMAX = W(1)
      IF (LBIG) BIG1 = THEMAX
C IDUP EQUALS ONE IF THERE WERE DUPLICATE ELEMENTS, ZERO OTHERWISE.
      IDUP = 0
      IF (IFLAG.EQ. (N+1)) IDUP = 1
      IF (IFLAG.LT.0) THEN
        IF (LP.NE.0) WRITE (LP,FMT=99994)
        GO TO 60

      END IF
C
C PERFORM ROW-GAUSS ELIMINATION ON THE STRUCTURE RECEIVED FROM MA28D/DD
      CALL MA30B(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IDISP,IKEEP(1,2),
     +           IKEEP(1,3),W,IW,IFLAG)
C
C TRANSFER COMMON BLOCK INFORMATION.
      IF (LBIG) BIG1 = BIG
      RMIN = MRMIN
      IF (IFLAG.GE.0) GO TO 50
      IFLAG = -2
      IF (LP.NE.0) WRITE (LP,FMT=99995)
      GO TO 60
C
C OPTIONALLY CALCULATE THE GROWTH PARAMETER.
   50 I1 = IDISP(1)
      IEND = LICN - I1 + 1
      IF (GROW) CALL MC24A(N,ICN,A(I1),IEND,IKEEP,IKEEP(1,4),W)
C INCREMENT ESTIMATE BY LARGEST ELEMENT IN INPUT MATRIX.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (GROW .AND. N.GT.1) W(2) = THEMAX
C SET FLAG IF THE ONLY ERROR IS DUE TO DUPLICATE ELEMENTS.
      IF (IDUP.EQ.1 .AND. IFLAG.GE.0) IFLAG = -14
      GO TO 70

   60 CONTINUE
   70 RETURN

99999 FORMAT (' ERROR RETURN FROM MA28B/BD WITH IFLAG=',I4,/,I7,
     +       ' ENTRIES DROPPED FROM STRUCTURE BY MA28A/AD')
99998 FORMAT (' ERROR RETURN FROM MA28B/BD BECAUSE ',
     +       'N OUT OF RANGE = ',I10)
99997 FORMAT (' ERROR RETURN FROM MA28B/BD BECAUSE ',
     +       'NZ NON POSITIVE = ',I10)
99996 FORMAT (' ERROR RETURN FROM MA28B/BD BECAUSE ',
     +       'LICN TOO SMALL = ',I10)
99995 FORMAT (' ERROR RETURN FROM MA28B/BD BECAUSE ',
     +       'ERROR RETURN FROM MA30B/BD')
99994 FORMAT (' ERROR RETURN FROM MA28B/BD ')

      END
      SUBROUTINE MA28C(N,A,LICN,ICN,IKEEP,RHS,W,MTYPE)
C
C THIS SUBROUTINE USES THE FACTORS FROM MA28A/AD OR MA28B/BD TO
C     SOLVE A SYSTEM OF EQUATIONS WITHOUT ITERATIVE REFINEMENT.
C THE PARAMETERS ARE ...
C N   INTEGER  ORDER OF MATRIX  NOT ALTERED BY SUBROUTINE.
C A   REAL/DOUBLE PRECISION ARRAY  LENGTH LICN.  THE SAME ARRAY AS
C     WAS USED IN THE MOST RECENT CALL TO MA28A/AD OR MA28B/BD.
C LICN  INTEGER  LENGTH OF ARRAYS A AND ICN.  NOT ALTERED BY
C     SUBROUTINE.
C ICN    INTEGER ARRAY OF LENGTH LICN.  SAME ARRAY AS OUTPUT FROM
C     MA28A/AD.  UNCHANGED BY MA28C/CD.
C IKEEP  INTEGER ARRAY OF LENGTH 5*N.  SAME ARRAY AS OUTPUT FROM
C     MA28A/AD.  UNCHANGED BY MA28C/CD.
C RHS    REAL/DOUBLE PRECISION ARRAY  LENGTH N.  ON ENTRY, IT HOLDS THE
C     RIGHT HAND SIDE.  ON EXIT, THE SOLUTION VECTOR.
C W      REAL/DOUBLE PRECISION ARRAY  LENGTH N. USED AS WORKSPACE BY
C     MA30C/CD.
C MTYPE  INTEGER  USED TO TELL MA30C/CD TO SOLVE THE DIRECT EQUATION
C     (MTYPE=1) OR ITS TRANSPOSE (MTYPE.NE.1).
C
C COMMON BLOCK VARIABLES.
C UNLESS OTHERWISE STATED COMMON BLOCK VARIABLES ARE AS IN MA28A/AD.
C     THOSE VARIABLES REFERENCED BY MA28C/CD ARE MENTIONED BELOW.
C RESID  REAL/DOUBLE PRECISION  VARIABLE RETURNS MAXIMUM RESIDUAL OF
C     EQUATIONS WHERE PIVOT WAS ZERO.
C MRESID  REAL/DOUBLE PRECISION VARIABLE USED BY MA28C/CD TO
C     COMMUNICATE BETWEEN MA28F/FD AND MA30H/HD.
C IDISP  INTEGER ARRAY  LENGTH 2  THE SAME AS THAT USED BY MA28A/AD.
C     IT IS UNCHANGED BY MA28B/BD.
C
C FURTHER INFORMATION ON COMMON BLOCK VARIABLES CAN BE FOUND IN BLOCK
C     DATA OR MA28A/AD.
C
C THIS CALL PERFORMS THE SOLUTION OF THE SET OF EQUATIONS.
C     .. Scalar Arguments ..
      INTEGER LICN,MTYPE,N
C     ..
C     .. Array Arguments ..
      REAL A(LICN),RHS(N),W(N)
      INTEGER ICN(LICN),IKEEP(N,5)
C     ..
C     .. External Subroutines ..
      EXTERNAL MA30C
C     ..
C     .. Common blocks ..
      COMMON /MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK,
     +       ABORT1,ABORT2
      COMMON /MA28G/IDISP(2)
      COMMON /MA30H/MRESID
      REAL EPS,MRESID,RESID,RMIN
      INTEGER ICNCP,IRANK,IRNCP,MINICN,MINIRN
      LOGICAL ABORT1,ABORT2
      INTEGER IDISP
C     ..
C     .. Save statement ..
      SAVE /MA28F/,/MA28G/,/MA30H/
C     ..
C     .. Executable Statements ..
      CALL MA30C(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IKEEP(1,5),IDISP,
     +           IKEEP(1,2),IKEEP(1,3),RHS,W,MTYPE)
C TRANSFER COMMON BLOCK INFORMATION.
      RESID = MRESID
      RETURN

      END
      SUBROUTINE MA28D(N,A,LICN,IVECT,JVECT,NZ,ICN,LENR,LENRL,LENOFF,IP,
     +                 IQ,IW1,IW,W1,IFLAG)
C THIS SUBROUTINE NEED NEVER BE CALLED BY THE USER DIRECTLY.
C     IT SORTS THE USER'S MATRIX INTO THE STRUCTURE OF THE DECOMPOSED
C     FORM AND CHECKS FOR THE PRESENCE OF DUPLICATE ENTRIES OR
C     NON-ZEROS LYING OUTSIDE THE SPARSITY PATTERN OF THE DECOMPOSITION
C     IT ALSO CALCULATES THE LARGEST ELEMENT IN THE INPUT MATRIX.
C     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0)
C     ..
C     .. Scalar Arguments ..
      REAL W1
      INTEGER IFLAG,LICN,N,NZ
C     ..
C     .. Array Arguments ..
      REAL A(LICN)
      INTEGER ICN(LICN),IP(N),IQ(N),IVECT(NZ),IW(N,2),IW1(N,3),
     +        JVECT(NZ),LENOFF(N),LENR(N),LENRL(N)
C     ..
C     .. Local Scalars ..
      REAL AA
      INTEGER I,IBLOCK,IDISP2,IDUMMY,II,INEW,IOLD,J1,J2,JCOMP,JDUMMY,JJ,
     +        JNEW,JOLD,MIDPT
      LOGICAL BLOCKL
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,AMAX1,IABS
C     ..
C     .. Common blocks ..
      COMMON /MA28E/LP,MP,LBLOCK,GROW
      COMMON /MA28G/IDISP(2)
      INTEGER LP,MP
      LOGICAL GROW,LBLOCK
      INTEGER IDISP
C     ..
C     .. Save statement ..
      SAVE /MA28E/,/MA28G/
C     ..
C     .. Executable Statements ..
      BLOCKL = LENOFF(1) .GE. 0
C IW1(I,3)  IS SET TO THE BLOCK IN WHICH ROW I LIES AND THE
C     INVERSE PERMUTATIONS TO IP AND IQ ARE SET IN IW1(.,1) AND
C     IW1(.,2) RESP.
C POINTERS TO BEGINNING OF THE PART OF ROW I IN DIAGONAL AND
C   OFF-DIAGONAL BLOCKS ARE SET IN IW(I,2) AND IW(I,1) RESP.
      IBLOCK = 1
      IW(1,1) = 1
      IW(1,2) = IDISP(1)
      DO 10 I = 1,N
        IW1(I,3) = IBLOCK
        IF (IP(I).LT.0) IBLOCK = IBLOCK + 1
        II = IABS(IP(I)+0)
        IW1(II,1) = I
        JJ = IQ(I)
        JJ = IABS(JJ)
        IW1(JJ,2) = I
        IF (I.EQ.1) GO TO 10
        IF (BLOCKL) IW(I,1) = IW(I-1,1) + LENOFF(I-1)
        IW(I,2) = IW(I-1,2) + LENR(I-1)
   10 CONTINUE
C PLACE EACH NON-ZERO IN TURN INTO ITS CORRECT LOCATION
C    IN THE A/ICN ARRAY.
      IDISP2 = IDISP(2)
      DO 170 I = 1,NZ
C NECESSARY TO AVOID REFERENCE TO UNASSIGNED ELEMENT OF ICN.
        IF (I.GT.IDISP2) GO TO 20
        IF (ICN(I).LT.0) GO TO 170
   20   IOLD = IVECT(I)
        JOLD = JVECT(I)
        AA = A(I)
C THIS IS A DUMMY LOOP FOR FOLLOWING A CHAIN OF INTERCHANGES.
C   IT WILL BE EXECUTED NZ TIMES IN TOTAL.
        DO 140 IDUMMY = 1,NZ
C PERFORM SOME VALIDITY CHECKS ON IOLD AND JOLD.
          IF (IOLD.LE.N .AND. IOLD.GT.0 .AND. JOLD.LE.N .AND.
     +        JOLD.GT.0) GO TO 30
          IF (LP.NE.0) WRITE (LP,FMT=99999) I,A(I),IOLD,JOLD
          IFLAG = -12
          GO TO 180

   30     INEW = IW1(IOLD,1)
          JNEW = IW1(JOLD,2)
C ARE WE IN A VALID BLOCK AND IS IT DIAGONAL OR OFF-DIAGONAL?
          IF (IW1(INEW,3)-IW1(JNEW,3)) 40,60,50
   40     IFLAG = -13
          IF (LP.NE.0) WRITE (LP,FMT=99998) IOLD,JOLD
          GO TO 180

   50     J1 = IW(INEW,1)
          J2 = J1 + LENOFF(INEW) - 1
          GO TO 110
C ELEMENT IS IN DIAGONAL BLOCK.
   60     J1 = IW(INEW,2)
          IF (INEW.GT.JNEW) GO TO 70
          J2 = J1 + LENR(INEW) - 1
          J1 = J1 + LENRL(INEW)
          GO TO 110

   70     J2 = J1 + LENRL(INEW)
C BINARY SEARCH OF ORDERED LIST  .. ELEMENT IN L PART OF ROW.
          DO 100 JDUMMY = 1,N
            MIDPT = (J1+J2)/2
            JCOMP = IABS(ICN(MIDPT)+0)
            IF (JNEW-JCOMP) 80,130,90
   80       J2 = MIDPT
            GO TO 100

   90       J1 = MIDPT
  100     CONTINUE
          IFLAG = -13
          IF (LP.NE.0) WRITE (LP,FMT=99997) IOLD,JOLD
          GO TO 180
C LINEAR SEARCH ... ELEMENT IN L PART OF ROW OR OFF-DIAGONAL BLOCKS.
  110     DO 120 MIDPT = J1,J2
            IF (IABS(ICN(MIDPT)+0).EQ.JNEW) GO TO 130
  120     CONTINUE
          IFLAG = -13
          IF (LP.NE.0) WRITE (LP,FMT=99997) IOLD,JOLD
          GO TO 180
C EQUIVALENT ELEMENT OF ICN IS IN POSITION MIDPT.
  130     IF (ICN(MIDPT).LT.0) GO TO 160
          IF (MIDPT.GT.NZ .OR. MIDPT.LE.I) GO TO 150
          W1 = A(MIDPT)
          A(MIDPT) = AA
          AA = W1
          IOLD = IVECT(MIDPT)
          JOLD = JVECT(MIDPT)
          ICN(MIDPT) = -ICN(MIDPT)
  140   CONTINUE
  150   A(MIDPT) = AA
        ICN(MIDPT) = -ICN(MIDPT)
        GO TO 170

  160   A(MIDPT) = A(MIDPT) + AA
C SET FLAG FOR DUPLICATE ELEMENTS.
        IFLAG = N + 1
  170 CONTINUE
C RESET ICN ARRAY  AND ZERO ELEMENTS IN L/U BUT NOT IN A.
C ALSO CALCULATE MAXIMUM ELEMENT OF A.
  180 W1 = ZERO
      DO 200 I = 1,IDISP2
        IF (ICN(I).LT.0) GO TO 190
        A(I) = ZERO
        GO TO 200

  190   ICN(I) = -ICN(I)
        W1 = AMAX1(W1,ABS(A(I)))
  200 CONTINUE
      RETURN

99999 FORMAT (' ELEMENT ',I6,' WITH VALUE ',1P,E12.4,/,10X,
     +       ' HAS INDICES ',I8,' ,',I8,' INDICES OUT OF RANGE')
99998 FORMAT (' NON-ZERO',I7,' ,',I6,' IN ZERO OFF-DIAGONAL BLOCK')
99997 FORMAT (' ELEMENT',I6,' ,',I6,' WAS NOT IN L/U PATTERN')

      END
      SUBROUTINE MA28I(N,NZ,AORG,IRNORG,ICNORG,LICN,A,ICN,IKEEP,RHS,X,R,
     +                 W,MTYPE,PREC,IFLAG)
C THIS SUBROUTINE USES THE FACTORS FROM AN EARLIER CALL TO MA28A/AD
C     OR MA28B/BD TO SOLVE THE SYSTEM OF EQUATIONS WITH ITERATIVE
C     REFINEMENT.
C
C THE PARAMETERS ARE...
C
C N IS EQUAL TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED BY THE
C     SUBROUTINE.
C NZ IS EQUAL TO THE NUMBER OF ENTRIES IN THE ORIGINAL MATRIX.  IT IS
C     NOT ALTERED BY THE SUBROUTINE.
C FOR THIS ENTRY THE ORIGINAL MATRIX MUST HAVE BEEN SAVED IN
C     AORG,IRNORG,ICNORG WHERE ENTRY AORG(K) IS IN ROW IRNORG(K) AND
C     COLUMN ICNORG(K), K=1,...NZ.  INFORMATION ABOUT THE FACTORS OF A
C     IS COMMUNICATED TO THIS SUBROUTINE VIA THE PARAMETERS LICN, A, ICN
C     AND IKEEP WHERE:
C AORG IS AN ARRAY OF LENGTH NZ.  NOT ALTERED BY MA28I/ID.
C IRNORG IS AN ARRAY OF LENGTH NZ.  NOT ALTERED BY MA28I/ID.
C ICNORG IS AN ARRAY OF LENGTH NZ.  NOT ALTERED BY MA28I/ID.
C LICN IS EQUAL TO THE LENGTH OF ARRAYS A AND ICN. IT IS NOT ALTERED BY
C     THE SUBROUTINE.
C A IS AN ARRAY OF LENGTH LICN. IT MUST BE UNCHANGED SINCE THE LAST CALL
C     TO MA28A/AD OR MA28B/BD. IT IS NOT ALTERED BY THE SUBROUTINE.
C ICN, IKEEP ARE THE ARRAYS (OF LENGTHS LICN AND 5*N, RESPECTIVELY) OF
C     THE SAME NAMES AS IN THE PREVIOUS ALL TO MA28A/AD. THEY SHOULD BE
C     UNCHANGED SINCE THIS EARLIER CALL AND    THEY ARE NOT ALTERED BY
C     MA28I/ID.
C THE  OTHER PARAMETERS ARE AS FOLLOWS:
C RHS IS AN ARRAY OF LENGTH N. THE USER MUST SET RHS(I) TO CONTAIN THE
C     VALUE OF THE I TH COMPONENT OF THE RIGHT HAND SIDE. IT IS NOT
C     ALTERED BY MA28I/ID.
C X IS AN ARRAY OF LENGTH N. IF AN INITIAL GUESS OF THE SOLUTION IS
C     GIVEN (ISTART EQUAL TO 1), THEN THE USER MUST SET X(I) TO CONTAIN
C     THE VALUE OF THE I TH COMPONENT OF THE ESTIMATED SOLUTION.  ON
C     EXIT, X(I) CONTAINS THE I TH COMPONENT OF THE SOLUTION VECTOR.
C R IS AN ARRAY OF LENGTH N. IT NEED NOT BE SET ON ENTRY.  ON EXIT, R(I)
C     CONTAINS THE I TH COMPONENT OF AN ESTIMATE OF THE ERROR IF MAXIT
C     IS GREATER THAN 0.
C W IS AN ARRAY OF LENGTH N. IT IS USED AS WORKSPACE BY MA28I/ID.
C MTYPE MUST BE SET TO DETERMINE WHETHER MA28I/ID WILL SOLVE A*X=RHS
C     (MTYPE EQUAL TO 1) OR AT*X=RHS (MTYPE NE 1, ZERO SAY). IT IS NOT
C     ALTERED BY MA28I/ID.
C PREC SHOULD BE SET BY THE USER TO THE RELATIVE ACCURACY REQUIRED. THE
C     ITERATIVE REFINEMENT WILL TERMINATE IF THE MAGNITUDE OF THE
C     LARGEST COMPONENT OF THE ESTIMATED ERROR RELATIVE TO THE LARGEST
C     COMPONENT IN THE SOLUTION IS LESS THAN PREC.  IT IS NOT ALTERED BY
C     MA28I/ID.
C IFLAG IS A DIAGNOSTIC FLAG WHICH WILL BE SET TO ZERO ON SUCCESSFUL
C     EXIT FROM MA28I/ID, OTHERWISE IT WILL HAVE A NON-ZERO VALUE. THE
C     NON-ZERO VALUE IFLAG CAN HAVE ON EXIT FROM MA28I/ID ARE ...
C     -16    INDICATING THAT MORE THAN MAXIT ITEARTIONS ARE REQUIRED.
C     -17    INDICATING THAT MORE CONVERGENCE WAS TOO SLOW.
C
C
C PRIVATE AND COMMON VARIABLES.
C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON.
C SEE COMMENTS IN CODE FOR USE OF PRIVATE VARIABLES.
C
C
C     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0)
C     ..
C     .. Scalar Arguments ..
      REAL PREC
      INTEGER IFLAG,LICN,MTYPE,N,NZ
C     ..
C     .. Array Arguments ..
      REAL A(LICN),AORG(NZ),R(N),RHS(N),W(N),X(N)
      INTEGER ICN(LICN),ICNORG(NZ),IKEEP(N,5),IRNORG(NZ)
C     ..
C     .. Local Scalars ..
      REAL CONVER,D,DD
      INTEGER I,ITERAT,NCOL,NROW
C     ..
C     .. External Subroutines ..
      EXTERNAL MA28C
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,AMAX1
C     ..
C     .. Common blocks ..
      COMMON /MA28E/LP,MP,LBLOCK,GROW
      COMMON /MA28H/TOL,THEMAX,BIG,DXMAX,ERRMAX,DRES,CGCE,NDROP,MAXIT,
     +       NOITER,NSRCH,ISTART,LBIG
      REAL BIG,CGCE,DRES,DXMAX,ERRMAX,THEMAX,TOL
      INTEGER ISTART,LP,MAXIT,MP,NDROP,NOITER,NSRCH
      LOGICAL GROW,LBIG,LBLOCK
C     ..
C     .. Save statement ..
      SAVE /MA28E/,/MA28H/
C     ..
C     .. Executable Statements ..
C
C   INITIALIZATION OF NOITER, ERRMAX AND IFLAG.
C
      NOITER = 0
      ERRMAX = ZERO
      DRES = ZERO
      IFLAG = 0
C
C  JUMP IF A STARTING VECTOR HAS BEEN SUPPLIED BY THE USER.
C
      IF (ISTART.EQ.1) GO TO 20
C
C  MAKE A COPY OF THE RIGHT-HAND SIDE VECTOR.
C
      DO 10 I = 1,N
        X(I) = RHS(I)
   10 CONTINUE
C
C  FIND THE FIRST SOLUTION.
C
      CALL MA28C(N,A,LICN,ICN,IKEEP,X,W,MTYPE)
C
C  STOP THE COMPUTATIONS IF   MAXIT=0.
C
   20 IF (MAXIT.EQ.0) GO TO 160
C
C  CALCULATE THE MAX-NORM OF THE FIRST SOLUTION.
C
      DD = 0.0
      DO 30 I = 1,N
        DD = AMAX1(DD,ABS(X(I)))
   30 CONTINUE
      DXMAX = DD
C
C  BEGIN THE ITERATIVE PROCESS.
C
      DO 120 ITERAT = 1,MAXIT
        D = DD
C
C  CALCULATE THE RESIDUAL VECTOR.
C
        DO 40 I = 1,N
          R(I) = RHS(I)
   40   CONTINUE
        IF (MTYPE.EQ.1) GO TO 60
        DO 50 I = 1,NZ
          NROW = IRNORG(I)
          NCOL = ICNORG(I)
          R(NCOL) = R(NCOL) - AORG(I)*X(NROW)
   50   CONTINUE
        GO TO 80
C MTYPE=1.
   60   DO 70 I = 1,NZ
          NROW = IRNORG(I)
          NCOL = ICNORG(I)
          R(NROW) = R(NROW) - AORG(I)*X(NCOL)
   70   CONTINUE
   80   DRES = 0.0
C
C  FIND THE MAX-NORM OF THE RESIDUAL VECTOR.
C
        DO 90 I = 1,N
          DRES = AMAX1(DRES,ABS(R(I)))
   90   CONTINUE
C
C  STOP THE CALCULATIONS IF THE MAX-NORM OF
C  THE RESIDUAL VECTOR IS ZERO.
C
        IF (DRES.EQ.0.0) GO TO 150
C
C  CALCULATE THE CORRECTION VECTOR.
C
        NOITER = NOITER + 1
        CALL MA28C(N,A,LICN,ICN,IKEEP,R,W,MTYPE)
C
C  FIND THE MAX-NORM OF THE CORRECTION VECTOR.
C
        DD = 0.0
        DO 100 I = 1,N
          DD = AMAX1(DD,ABS(R(I)))
  100   CONTINUE
C
C  CHECK THE CONVERGENCE.
C
        IF (DD.GT.D*CGCE .AND. ITERAT.GE.2) GO TO 130
        IF (DXMAX*10.0+DD.EQ.DXMAX*10.0) GO TO 140
C
C  ATTEMPT TO IMPROVE THE SOLUTION.
C
        DXMAX = 0.0
        DO 110 I = 1,N
          X(I) = X(I) + R(I)
          DXMAX = AMAX1(DXMAX,ABS(X(I)))
  110   CONTINUE
C
C  CHECK THE STOPPING CRITERION.
C
        IF (DD.LT.PREC*DXMAX) GO TO 140
  120 CONTINUE
C MORE THAN MAXIT ITERATIONS REQUIRED.
      IFLAG = -16
      WRITE (LP,FMT=99999) IFLAG,MAXIT
      GO TO 140
C CONVERGENCE RATE UNACCEPTABLY SLOW.
  130 IFLAG = -17
      CONVER = DD/D
      WRITE (LP,FMT=99998) IFLAG,CONVER,CGCE
C
C  THE ITERATIVE PROCESS IS TERMINATED.
C
  140 ERRMAX = DD
  150 CONTINUE
  160 RETURN

99999 FORMAT (' ERROR RETURN FROM MA28I/ID WITH IFLAG = ',I3,/,' MORE ',
     +       'THAN',I5,' ITERATIONS REQUIRED')
99998 FORMAT (' ERROR RETURN FROM MA28I WITH IFLAG = ',I3,/,' CONVERGE',
     +       'NCE RATE OF ',1P,E9.2,' TOO SLOW',/,
     +       ' MAXIMUM ACCEPTABLE RATE',' SET TO ',1P,E9.2)

      END
      BLOCK DATA MA28J
C
C COMMENTS ON ALL THE COMMON BLOCK VARIABLES ARE GIVEN HERE EVEN
C     THOUGH SOME ARE NOT INITIALIZED BY BLOCK DATA.
C LP,MP ARE USED BY THE SUBROUTINE AS THE UNIT NUMBERS FOR ITS WARNING
C     AND DIAGNOSTIC MESSAGES. DEFAULT VALUE FOR BOTH IS 6 (FOR LINE
C     PRINTER OUTPUT). THE USER CAN EITHER RESET THEM TO A DIFFERENT
C     STREAM NUMBER OR SUPPRESS THE OUTPUT BY SETTING THEM TO ZERO.
C     WHILE LP DIRECTS THE OUTPUT OF ERROR DIAGNOSTICS FROM THE
C     PRINCIPAL SUBROUTINES AND INTERNALLY CALLED SUBROUTINES, MP
C     CONTROLS ONLY THE OUTPUT OF A MESSAGE WHICH WARNS THE USER THAT HE
C     HAS INPUT TWO OR MORE NON-ZEROS A(I), . . ,A(K) WITH THE SAME ROW
C     AND COLUMN INDICES.  THE ACTION TAKEN IN THIS CASE IS TO PROCEED
C     USING A NUMERICAL VALUE OF A(I)+...+A(K). IN THE ABSENCE OF OTHER
C     ERRORS, IFLAG WILL EQUAL -14 ON EXIT.
C LBLOCK IS A LOGICAL VARIABLE WHICH CONTROLS AN OPTION OF FIRST
C     PREORDERING THE MATRIX TO BLOCK LOWER TRIANGULAR FORM (USING
C     HARWELL SUBROUTINE MC23A). THE PREORDERING IS PERFORMED IF LBLOCK
C     IS EQUAL TO ITS DEFAULT VALUE OF .TRUE. IF LBLOCK IS SET TO
C     .FALSE. , THE OPTION IS NOT INVOKED AND THE SPACE ALLOCATED TO
C     IKEEP CAN BE REDUCED TO 4*N+1.
C GROW IS A LOGICAL VARIABLE. IF IT IS LEFT AT ITS DEFAULT VALUE OF
C     .TRUE. , THEN ON RETURN FROM MA28A/AD OR MA28B/BD, W(1) WILL GIVE
C     AN ESTIMATE (AN UPPER BOUND) OF THE INCREASE IN SIZE OF ELEMENTS
C     ENCOUNTERED DURING THE DECOMPOSITION. IF THE MATRIX IS WELL
C     SCALED, THEN A HIGH VALUE FOR W(1), RELATIVE TO THE LARGEST ENTRY
C     IN THE INPUT MATRIX, INDICATES THAT THE LU DECOMPOSITION MAY BE
C     INACCURATE AND THE USER SHOULD BE WARY OF HIS RESULTS AND PERHAPS
C     INCREASE U FOR SUBSEQUENT RUNS.  WE WOULD LIKE TO EMPHASISE THAT
C     THIS VALUE ONLY RELATES TO THE ACCURACY OF OUR LU DECOMPOSITION
C     AND GIVES NO INDICATION AS TO THE SINGULARITY OF THE MATRIX OR THE
C     ACCURACY OF THE SOLUTION.  THIS UPPER BOUND CAN BE A SIGNIFICANT
C     OVERESTIMATE PARTICULARLY IF THE MATRIX IS BADLY SCALED. IF AN
C     ACCURATE VALUE FOR THE GROWTH IS REQUIRED, LBIG (Q.V.) SHOULD BE
C     SET TO .TRUE.
C EPS,RMIN ARE REAL VARIABLES. IF, ON ENTRY TO MA28B/BD, EPS IS LESS
C     THAN ONE, THEN RMIN WILL GIVE THE SMALLEST RATIO OF THE PIVOT TO
C     THE LARGEST ELEMENT IN THE CORRESPONDING ROW OF THE UPPER
C     TRIANGULAR FACTOR THUS MONITORING THE STABILITY OF SUCCESSIVE
C     FACTORIZATIONS. IF RMIN BECOMES VERY LARGE AND W(1) FROM
C     MA28B/BD IS ALSO VERY LARGE, IT MAY BE ADVISABLE TO PERFORM A
C     NEW DECOMPOSITION USING MA28A/AD.
C RESID IS A REAL VARIABLE WHICH ON EXIT FROM MA28C/CD GIVES THE VALUE
C     OF THE MAXIMUM RESIDUAL OVER ALL THE EQUATIONS UNSATISFIED BECAUSE
C     OF DEPENDENCY (ZERO PIVOTS).
C IRNCP,ICNCP ARE INTEGER VARIABLES WHICH MONITOR THE ADEQUACY OF "ELBOW
C     ROOM" IN IRN AND A/ICN RESPECTIVELY. IF EITHER IS QUITE LARGE (SAY
C     GREATER THAN N/10), IT WILL PROBABLY PAY TO INCREASE THE SIZE OF
C     THE CORRESPONDING ARRAY FOR SUBSEQUENT RUNS. IF EITHER IS VERY LOW
C     OR ZERO THEN ONE CAN PERHAPS SAVE STORAGE BY REDUCING THE SIZE OF
C     THE CORRESPONDING ARRAY.
C MINIRN,MINICN ARE INTEGER VARIABLES WHICH, IN THE EVENT OF A
C     SUCCESSFUL RETURN (IFLAG GE 0 OR IFLAG=-14) GIVE THE MINIMUM SIZE
C     OF IRN AND A/ICN RESPECTIVELY WHICH WOULD ENABLE A SUCCESSFUL RUN
C     ON AN IDENTICAL MATRIX. ON AN EXIT WITH IFLAG EQUAL TO -5, MINICN
C     GIVES THE MINIMUM VALUE OF ICN FOR SUCCESS ON SUBSEQUENT RUNS ON
C     AN IDENTICAL MATRIX. IN THE EVENT OF FAILURE WITH IFLAG= -6, -4,
C     -3, -2, OR -1, THEN MINICN AND MINIRN GIVE THE MINIMUM VALUE OF
C     LICN AND LIRN RESPECTIVELY WHICH WOULD BE REQUIRED FOR A
C     SUCCESSFUL DECOMPOSITION UP TO THE POINT AT WHICH THE FAILURE
C     OCCURRED.
C IRANK IS AN INTEGER VARIABLE WHICH GIVES AN UPPER BOUND ON THE RANK OF
C     THE MATRIX.
C ABORT1 IS A LOGICAL VARIABLE WITH DEFAULT VALUE .TRUE.  IF ABORT1 IS
C     SET TO .FALSE.  THEN MA28A/AD WILL DECOMPOSE STRUCTURALLY SINGULAR
C     MATRICES (INCLUDING RECTANGULAR ONES).
C ABORT2 IS A LOGICAL VARIABLE WITH DEFAULT VALUE .TRUE.  IF ABORT2 IS
C     SET TO .FALSE. THEN MA28A/AD WILL DECOMPOSE NUMERICALLY SINGULAR
C     MATRICES.
C IDISP IS AN INTEGER ARRAY OF LENGTH 2. ON OUTPUT FROM MA28A/AD, THE
C     INDICES OF THE DIAGONAL BLOCKS OF THE FACTORS LIE IN POSITIONS
C     IDISP(1) TO IDISP(2) OF A/ICN. THIS ARRAY MUST BE PRESERVED
C     BETWEEN A CALL TO MA28A/AD AND SUBSEQUENT CALLS TO MA28B/BD,
C     MA28C/CD OR MA28I/ID.
C TOL IS A REAL VARIABLE.  IF IT IS SET TO A POSITIVE VALUE, THEN ANY
C     NON-ZERO WHOSE MODULUS IS LESS THAN TOL WILL BE DROPPED FROM THE
C     FACTORIZATION.  THE FACTORIZATION WILL THEN REQUIRE LESS STORAGE
C     BUT WILL BE INACCURATE.  AFTER A RUN OF MA28A/AD WITH TOL POSITIVE
C     IT IS NOT POSSIBLE TO USE MA28B/BD AND THE USER IS RECOMMENDED TO
C     USE MA28I/ID TO OBTAIN THE SOLUTION.  THE DEFAULT VALUE FOR TOL IS
C     0.0.
C THEMAX IS A REAL VARIABLE.  ON EXIT FROM MA28A/AD, IT WILL HOLD THE
C     LARGEST ENTRY OF THE ORIGINAL MATRIX.
C BIG IS A REAL VARIABLE. IF LBIG HAS BEEN SET TO .TRUE., BIG WILL HOLD
C     THE LARGEST ENTRY ENCOUNTERED DURING THE FACTORIZATION BY MA28A/AD
C     OR MA28B/BD.
C DXMAX IS A REAL VARIABLE. ON EXIT FROM MA28I/ID, DXMAX WILL BE SET TO
C     THE LARGEST COMPONENT OF THE SOLUTION.
C ERRMAX IS A REAL VARIABLE.  ON EXIT FROM MA28I/ID, IF MAXIT IS
C     POSITIVE, ERRMAX WILL BE SET TO THE LARGEST COMPONENT IN THE
C     ESTIMATE OF THE ERROR.
C DRES IS A REAL VARIABLE.  ON EXIT FROM MA28I/ID, IF MAXIT IS POSITIVE,
C     DRES WILL BE SET TO THE LARGEST COMPONENT OF THE RESIDUAL.
C CGCE IS A REAL VARIABLE. IT IS USED BY MA28I/ID TO CHECK THE
C     CONVERGENCE RATE.  IF THE RATIO OF SUCCESSIVE CORRECTIONS IS
C     NOT LESS THAN CGCE THEN WE TERMINATE SINCE THE CONVERGENCE
C     RATE IS ADJUDGED TOO SLOW.
C NDROP IS AN INTEGER VARIABLE. IF TOL HAS BEEN SET POSITIVE, ON EXIT
C     FROM MA28A/AD, NDROP WILL HOLD THE NUMBER OF ENTRIES DROPPED FROM
C     THE DATA STRUCTURE.
C MAXIT IS AN INTEGER VARIABLE. IT IS THE MAXIMUM NUMBER OF ITERATIONS
C     PERFORMED BY MA28I/ID. IT HAS A DEFAULT VALUE OF 16.
C NOITER IS AN INTEGER VARIABLE. IT IS SET BY MA28I/ID TO THE NUMBER OF
C     ITERATIVE REFINEMENT ITERATIONS ACTUALLY USED.
C NSRCH IS AN INTEGER VARIABLE. IF NSRCH IS SET TO A VALUE LESS THAN N,
C     THEN A DIFFERENT PIVOT OPTION WILL BE EMPLOYED BY MA28A/AD.  THIS
C     MAY RESULT IN DIFFERENT FILL-IN AND EXECUTION TIME FOR MA28A/AD.
C     IF NSRCH IS LESS THAN OR EQUAL TO N, THE WORKSPACE ARRAY IW CAN BE
C     REDUCED IN LENGTH.  THE DEFAULT VALUE FOR NSRCH IS 32768.
C ISTART IS AN INTEGER VARIABLE. IF ISTART IS SET TO A VALUE OTHER THAN
C     ZERO, THEN THE USER MUST SUPPLY AN ESTIMATE OF THE SOLUTION TO
C     MA28I/ID.  THE DEFAULT VALUE FOR ISTART IS ZERO.
C LBIG IS A LOGICAL VARIABLE. IF LBIG IS SET TO .TRUE., THE VALUE OF THE
C     LARGEST ELEMENT ENCOUNTERED IN THE FACTORIZATION BY MA28A/AD OR
C     MA28B/BD IS RETURNED IN BIG.  SETTING LBIG TO .TRUE.  WILL
C     INCREASE THE TIME FOR MA28A/AD MARGINALLY AND THAT FOR MA28B/BD
C     BY ABOUT 20%.  THE DEFAULT VALUE FOR LBIG IS .FALSE.
C
C     COMMON /MA28G/ IDISP(2)
C     .. Common blocks ..
      COMMON /MA28E/LP,MP,LBLOCK,GROW
      COMMON /MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK,
     +       ABORT1,ABORT2
      COMMON /MA28H/TOL,THEMAX,BIG,DXMAX,ERRMAX,DRES,CGCE,NDROP,MAXIT,
     +       NOITER,NSRCH,ISTART,LBIG
      REAL BIG,CGCE,DRES,DXMAX,EPS,ERRMAX,RESID,RMIN,THEMAX,TOL
      INTEGER ICNCP,IRANK,IRNCP,ISTART,LP,MAXIT,MINICN,MINIRN,MP,NDROP,
     +        NOITER,NSRCH
      LOGICAL ABORT1,ABORT2,GROW,LBIG,LBLOCK
C     ..
C     .. Save statement ..
      SAVE /MA28E/,/MA28F/,/MA28H/
C     ..
C     .. Data statements ..
      DATA EPS/1.0E-4/,TOL/0.0E0/,CGCE/0.5E0/
      DATA MAXIT/16/
      DATA LP/6/,MP/6/,NSRCH/32768/,ISTART/0/
      DATA LBLOCK/.TRUE./,GROW/.TRUE./,LBIG/.FALSE./
      DATA ABORT1/.TRUE./,ABORT2/.TRUE./
C     ..
C     .. Executable Statements ..
      END
